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ABSTRACT 


This  paper  discusses  the  modeling  of  non-equilibrium  ef¬ 
fects  in  inductively  coupled  plasma  facilities.  The  model 
relies  on  the  solution  of  the  Navier-Stokes  and  Maxwell 
equations  in  a  one-dimensional  geometry.  Steady-state 
solutions  are  obtained  by  means  of  an  implicit  Finite 
Volume  method.  Non-equilibrium  effects  are  treated  by 
means  of  a  hybrid  State-to-State  formulation.  The  elec¬ 
tronic  states  of  atoms  are  treated  as  separate  species,  al¬ 
lowing  for  non-Boltzmann  distributions  of  their  popula¬ 
tions.  Thermal  non-equilibrium  between  the  translation 
and  vibrational  of  heavy-particles  is  accounted  for  by 
means  of  a  multi-temperature  approach.  The  results  show 
that  non-equilibrium  plays  an  important  role  close  to  the 
walls,  due  to  the  combined  effects  of  Ohming  heating, 
and  chemical  composition  and  temperature  gradients. 

Key  words:  NLTE  plasmas,  State-to-State  modeling.  In¬ 
ductively  Coupled  Plasma  torch. 


1.  INTRODUCTION 

Inductively  coupled  plasma  (ICP)  torches  have  wide 
range  of  possible  applications  which  include  deposition 
of  metal  coatings,  synthesis  of  ultra-fine  powders,  gener¬ 
ation  of  high  purity  silicon  and  testing  of  thermal  protec¬ 
tion  materials  for  atmospheric  entry  vehicles  [1,2]. 

In  its  simplest  configuration,  an  ICP  torch  consists  of 
a  quartz  tube  surrounded  by  an  inductor  coil  made  of 
a  series  of  parallel  current-carrying  rings.  The  radio¬ 
frequency  (RF)  currents  running  through  the  inductor  in¬ 
duce  toroidal  currents  in  the  gas  which  is  heated  thanks  to 
Ohmic  dissipation  [2,  3],  If  the  energy  supplied  is  large 
enough,  the  gas  flowing  through  the  torch  can  undergo 
ionization,  leading  to  the  formation  of  a  plasma. 

The  accurate  modeling  of  the  flow-field  and  electromag¬ 
netic  phenomena  inside  an  ICP  torch  requires  the  cou¬ 
pled  solution  of  the  Navier-Stokes  and  Maxwell  equa¬ 
tions  [2].  In  literature,  the  Local  Thermodynamic  Equi¬ 
librium  (LTE)  assumption  is  often  used  to  describe  the 
state  of  the  gas  in  the  discharge  region  [4-17],  However, 


Non  Local  Thermodynamic  Equilibrium  (NLTE)  simula¬ 
tions  of  Argon  [18,  19]  and  air  plasmas  [20],  have  shown 
that  the  LTE  assumption  may  not  always  hold. 

An  accurate  modeling  of  NLTE  effects  in  ICP  RF  plasmas 
can  be  achieved  by  means  of  State-to-State  models  [21- 
33].  These  treat  each  internal  energy  state  as  a  separate 
pseudo  species,  thus  allowing  for  non-Boltzmann  distri¬ 
butions.  Rate  coefficients  are  usually  obtained  through 
quantum  chemistry  calculations  [34-39]  or  through  phe¬ 
nomenological  models  [40,  41].  State-to-State  models 
provide  a  superior  description  compared  to  conventional 
multi-temperature  models,  which  are  based  on  Maxwell- 
Boltzmann  distributions  [  42 — 45  ] .  However,  due  to  the 
large  number  of  governing  equations  to  be  solved,  their 
application  to  multi-dimensional  problems  can  become 
computationally  demanding  [46-50]. 

The  purpose  of  the  present  work  is  the  development  of 
a  magneto-hydrodynamic  NLTE  model  for  an  ICP  torch. 
The  flow  model  has  to  be  simple  enough  to  allow  for  the 
implementation  of  sophisticated  NLTE  models  describ¬ 
ing  the  kinetics  of  atoms  and  molecules  in  the  inductor 
region.  The  final  goal  is  the  assessment  of  the  extent  of 
non-equilibrium  phenomena  occurring  in  ICP  generators 
used  for  Thermal  Protection  System  (TPS)  testing. 

The  paper  is  structured  as  follows.  Section  2  describes 
the  physical  model.  The  numerical  method  is  presented 
in  Sect.  3.  Computational  results  are  discussed  in  Sect. 
4.  Conclusions  are  outlined  in  Sect.  5. 

2.  PHYSICAL  MODELING 

The  NLTE  model  for  the  ICP  torch  is  built  based  on  the 
torch  geometry  displayed  in  Figure  1 .  To  make  the  prob¬ 
lem  tractable,  the  following  assumptions  are  introduced: 

(i)  Constant  pressure  and  no  macroscopic  streaming, 

(ii)  Charge  neutrality  and  no  displacement  current, 

(iii)  Steady-state  conditions  for  gas  quantities  (i.e., 
d()/dt  =  0), 

(iv)  No  gradients  along  the  axial  and  circumferential  di¬ 
rections  (i.e.,  d()/dz  =  0,  d()/d(f>  =  0). 
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Figure  1:  Torch  geometry  and  adopted  reference  frame. 


2.1.  Electromagnetic  field 


The  electromagnetic  field  inside  the  ICP  torch  is  de¬ 
scribed  by  the  Maxwell  equations: 

V  ■  E  =  —  (1) 

eo 

V  •  B  =  0  (2) 

9B 

VxE=-—  (3) 

at 

<9E 

V  x  B  =  p03  +  Moeo^r  (4) 

at 

where  quantities  E  and  B  are  the  electric  and  magnetic 
fields,  respectively.  Quantity  pc  stands  for  the  charge  den¬ 
sity.  The  current  density  J  is  assumed  to  obey  Ohm’s 
law  J  =  treE  [51],  with  ae  being  the  electrical  conduc¬ 
tivity.  Quantities  e0  and  p0  are  the  vacuum  permittiv¬ 
ity  and  magnetic  permeability,  respectively.  The  applica¬ 
tion  of  the  simplifying  assumptions  just  introduced  to  the 
Maxwell  equations  (l)-(4)  leads  to  the  following  induc¬ 
tion  equation  for  the  induced  toroidal  electric  field: 


d  (l  drE^  \  _  dE^ 
dr  \r  dr  )  ~  M°CTe  dt 


(5) 


Since  the  induced  eddy  currents  which  are  responsible  for 
the  heating  of  the  gas  are  induced  by  a  primary  current 
whose  intensity  varies  sinusoidally  in  time,  it  seems  nat¬ 
ural  to  seek  for  a  mono-chromatic  wave  solution,  E,p  = 
E  exp(*27r/f),  where  /  is  the  frequency  of  the  primary 
current.  To  account  for  the  possible  phase  difference  be¬ 
tween  the  induced  electric  and  magnetic  fields,  the  am¬ 
plitude  E  is  taken  complex,  E  =  Ere  +  iEun.  The  substi¬ 
tution  of  E  exp(%2ir ft)  in  Eq.  (5)  leads  to: 
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The  electromagnetic  (em)  conservative  variable,  flux  and 
source  term  vectors  are: 


Uem  =  [  E„ 
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where  the  angular  frequency  is  u  =  2tt/.  The  symbol  T 
denotes  the  transpose  operator. 


Equation  (6)  must  be  supplemented  with  boundary  condi¬ 
tions  at  the  axis  (r  =  0)  and  at  the  torch  wall  (r  =  R,  with 
R  being  the  torch  radius).  On  the  axis,  due  to  symmetry, 
both  the  real  and  imaginary  components  of  the  induced 
electric  field  must  vanish: 

Eie  =  0,  Eim  =  0,  at  r  =  0  (10) 


The  boundary  condition  at  the  torch  wall  is  obtained  as 
follows.  The  amplitudes  of  the  toroidal  electric  field  and 
the  axial  magnetic  field  are  linked  via  [17]: 


1  drE 
r  dr 


=  —uoB 


(ID 


where  the  amplitude  B  is  taken  complex.  Immediately 
outside  the  wall  the  magnetic  field  must  be  real  and,  since 
there  is  no  plasma  outside  the  tube,  its  value  can  only 
depend  on  the  ICP  operating  conditions  and  characteris¬ 
tics.  If  the  torch  is  long  enough,  the  magnetic  field  at  the 
torch  wall  can  be  approximated  with  the  expression  for 
an  infinite  solenoid,  B{r  =  R)  =  po NIC,  where  quan¬ 
tities  N  and  Ic  are  the  number  of  turns  per  unit-length 
and  the  amplitude  of  the  primary  current.  The  evaluation 
of  Eq.  ( 1 1 )  at  the  torch  wall  and  the  use  of  the  relation 
B(r  =  R)  =  PqNIc  gives  the  wall  boundary  condition 
for  the  induced  electric  field: 


1  drEle  _  o  1  dr  Ejm 

r  dr  ’  r  dr 


— u>pqNIc ,  at  r  =  R 

(12) 


2.2.  Gas  field 

The  gas  contained  in  the  torch  is  made  of  electrons, 
atoms  and  (diatomic)  molecules.  Charged  particles  com¬ 
prise  electrons  and  positively  singly  charged  ions.  The 
set  S  stores  the  chemical  components,  and  the  heavy- 
particles  are  stored  in  the  set  5),  ■  The  atomic  and  molec¬ 
ular  components  are  stored  in  the  sets  S-d  and  iSm,  respec¬ 
tively.  The  previously  introduced  sets  satisfy  the  relations 
<Sh  =  S-d  U  5m  and  S  =  {e~}  U  <Sm,  where  the  symbol  e- 
indicates  the  free-electrons.  The  electronic  levels  of  the 
heavy  components  are  stored  in  sets  if  (with  s  £  iSh )  and 
are  treated  as  separate  pseudo-species  based  on  a  State- 
to-State  approach  [52].  The  notation  s,:  is  then  introduced 
to  denote  the  i-th  electronic  level  of  the  heavy  component 
s  £  <Sh,  with  the  related  degeneracy  and  energy  being 
gf.  and  Ef  ,  respectively.  A  multi-temperature  model  is 
instead  used  for  vibration  of  molecules  and  translation 
of  free-electrons  (with  the  related  temperatures  being  Tv 
and  Te,  respectively)  [53].  Rotational  non-equilibrium 
effects  are  disregarded. 


Thermodynamics 

The  gas  pressure  is  computed  based  on  Dalton’s  law  of 
partial  pressures  by  summing  the  contributions  of  free- 
electrons  and  heavy-particles,  p  =  pe  +  pn,  where  the 


symbol  k|>  stands  for  Boltzmann’s  constant.  The  partial 
pressures  of  free-electrons  and  heavy-particles  are,  re¬ 
spectively,  pe  =  nek^Te  and  p h  =  n^k^T,  where  quan¬ 
tities  ne  and  rih  denote,  respectively,  the  related  the  num¬ 
ber  densities.  The  number  density  of  heavy-particles  is 
obtained  via  nh  =  Y2s&sh  n°  >  with  ns  =  Dieij1  n 

The  gas  total,  vibrational  and  free-electron  energy  densi¬ 
ties  read: 

Pe  =  \v+  5>[2SCT)  +  K{Tw)}  +  J2nsAEl+ 

sG«Sh 

E  C3) 

sG«Sh 

iexf 

pev  =  Ens^(Tv),  pee  =  (14) 

Quantity  A Els  stands  for  the  formation  energy  (per  parti¬ 
cle)  of  the  heavy  component  s  £  iSh-  The  average  particle 
rotational  and  vibrational  energies  (Ers  and  Eg,  s  £  Sm 
respectively)  are  computed,  respectively,  according  to 
the  rigid-rotor  and  harmonic-oscillator  models  [54] . 

Thermodynamic  data  used  in  this  work  are  taken  from 
Gurvich  tables  [55]  (with  the  exception  of  the  spectro¬ 
scopic  data  for  the  electronic  levels  taken  from  [26]). 


Kinetics 

Collisional  processes  The  NLTE  model  for  ICP  RF 
plasmas  developed  in  this  work  accounts  for  the  follow¬ 
ing  kinetic  processes: 

(i)  Excitation  by  electron  impact  (EXe), 

(ii)  Ionization  by  electron  impact  (Ie), 

(iii)  Dissociation  by  electron  impact  (De), 

(iv)  Dissociation  by  heavy-particle  impact  (Dh), 

(v)  Associative  ionization  (AI). 


Rate  coefficients  and  production  terms  The  en¬ 
dothermic  rate  coefficients  for  electron  induced  processes 
and  associative  ionization  reactions  are  taken  from  the 
ABBA  model  [26-29].  Those  for  dissociation  by  heavy- 
particle  impact  are  taken  from  the  work  of  Park  [43]. 
Reverse  rate  coefficients  are  obtained  based  on  micro¬ 
reversibility  [56,  57]. 

The  mass  production  terms  for  free-electrons  and  heavy- 
particles  due  to  the  kinetic  processes  considered  in  this 
work  are  computed  based  on  the  zeroth-order  reaction 
rate  theory  [56,  57].  In  what  follows,  the  latter  quanti¬ 
ties  are  indicated  with  the  notation  oje  and  uiSi . 


The  energy  transfer  terms  for  the  gas  vibrational  en¬ 
ergy  account  for  (i)  vibrational-translational  (vt)  en¬ 
ergy  exchange  in  molecule  heavy-particle  collisions,  (ii) 
vibrational-electron  (ve)  energy  exchange  in  molecule 
electron  collisions,  and  (iii)  the  creation/destruction  of  vi¬ 
brational  energy  in  chemical  reactions  (cv).  The  first  two 
energy  transfer  terms  (indicated  in  what  follows  with  1 2vt 
and  Ove,  respectively)  are  evaluated  based  on  a  Landau- 
Teller  model  [58],  while  the  chemistry- vibration  coupling 
term  (Ocv)  is  computed  by  using  the  non-preferential  dis¬ 
sociation  model  of  Candler  [59].  The  relaxation  times  for 
vt  energy  transfer  are  computed  by  means  of  the  mod¬ 
ified  formula  of  Millikan  and  White  proposed  by  Park 
[43].  The  energy  transfer  in  molecule-electron  inelastic 
collisions  is  considered  only  for  N2.  The  corresponding 
relaxation  time  is  taken  from  the  work  of  Bourdon  [60]. 

The  energy  transfer  terms  for  the  free-electron  gas  ac¬ 
count  for  energy  exchange  undergone  by  free-electrons 
in  (i)  elastic  collisions  with  heavy-particles  (Dei),  (ii)  in¬ 
elastic  electron  induced  excitation,  ionization  and  disso¬ 
ciation  processes  (Di„)  and  (iii)  Joule  heating  (Dj).  The 
expressions  for  the  first  two  can  be  found  in  [27-29].  The 
(time-averaged)  Joule  heating  source  term  is  obtained  by 
averaging  over  a  period  the  instantaneous  Joule  heating 
power  and  reads  [15,  17]: 

^=l-ae{El  +  El >m)  (15) 


Transport 

Transport  phenomena  are  treated  based  on  the  results 
obtained  from  the  application  of  the  Chapman-Enskog 
method  to  the  Boltzmann  equation  [61].  In  the  present 
work,  transport  phenomena  are  modeled  by  assuming  that 
(i)  inelastic  and  reactive  collisions  have  a  no  effect  on  the 
transport  properties  and  fluxes  and  (ii)  the  collision  cross- 
sections  for  elastic  scattering  do  not  depend  on  the  inter¬ 
nal  quantum  states.  In  view  of  the  last  assumption,  the  di¬ 
mension  of  the  transport  matrices  shrinks  from  the  num¬ 
ber  of  species  to  the  number  of  chemical  components. 


Heavy-particle  gas  The  translational  component  of 
thermal  conductivity  (At)  is: 

At  =  E a*x>  (16) 

s(E«Sh 

where  the  mole  fractions  of  the  heavy  components  are 
Xs  =  risksT Ip  ( s  6  <Sh).  The  coefficients  a*  are  solu¬ 
tion  of  the  linear  (symmetric)  transport  system: 

E  G =  Xs  (17) 

p&Sh 

s  €  <Sh-  Quantities  G^p  are  the  entries  of  the  thermal 
conductivity  (symmetric)  transport  matrix  [56].  The  con¬ 
tributions  of  the  gas  rotational  and  vibrational  degrees  of 


freedom  to  the  thermal  conductivity  (Ar  and  Av,  respec¬ 
tively)  are  taken  into  account  by  means  of  the  generalized 
Eucken’s  correction  [56]. 


Electron  gas  The  thermal  and  electrical  conductivity 
of  the  electron  gas  are  (in  the  second  and  third-order 
Laguerre-Sonine  approximations,  respectively)  [62,  63]: 


/27rkBTe  XeA2J 

(18) 

/  me  A£A“-(A£)2 

3  e2 

/  27rkB  XeA 11 

(19) 

2  kB  V 

meTe  A°e°A ii  -  (A ie0)2 

where  the  mole  fraction  of  free-electrons  is  Xe  = 
nekBTe/p  and  e  =  1.602  x  10~19C  is  the  electron 
charge.  Quantities  Ale  are  the  Devoto  collision  integrals 
[62], 


where  the  gas  (g)  conservative  variable,  flux  and  source 
term  vectors  are: 


Ug  = 

[  Pe 

PSi 

pe 

„  ,T 

p6y  pec  J 

(24) 
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[  Je 

Js, 
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qe  \ 

(25) 

Sg  = 

[  We 

wSi 

Oj 

fly  fle  ]T 

(26) 

i  £  If1,  s  £  <Sh,  with  flv  =  flvt  +  f2ve  +  flcv  and  fle  = 
flel  +  flin  T  • 

The  boundary  conditions  used  for  solving  Eq.  (23)  are  a 
symmetry  boundary  condition  at  the  axis  and  an  isother¬ 
mal  non-catalytic  boundary  condition  at  the  torch  wall 
(with  the  related  imposed  wall  temperature  indicated  with 
Tw).  In  analogy  with  the  work  of  Mostaghimi  et  al. 
[18,  19],  an  adiabatic  boundary  condition  is  used  for  the 
vibrational  and  free-electron  temperatures  at  the  torch 
wall. 


3.  NUMERICAL  METHOD 


Mass  diffusion  and  heat  fluxes  The  mass  diffusion 
fluxes  are  found  by  solving  the  Stefan-Maxwell  equa¬ 
tions  under  the  constraints  of  global  mass  conservation 
and  ambipolar  diffusion  [62-65].  The  diffusion  driving 
forces  include  only  mole  fraction  gradients  (thermal  and 
baro  diffusion  are  both  neglected).  In  view  of  the  as¬ 
sumed  independence  of  the  elastic  collision  cross-section 
on  the  internal  quantum  states,  the  Stefan-Maxwell  equa¬ 
tions  are  solved  for  the  free-electron  and  heavy  compo¬ 
nent  diffusion  fluxes  (Je  and  Js,  s  £  <Si,,  respectively). 
The  mass  diffusion  fluxes  for  the  internal  (electronic)  lev¬ 
els  are  then  found  as  shown  in  [17].  The  total,  vibrational 
and  free-electron  heat  fluxes  are: 
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3.1.  Coupled  formulation 


The  governing  equations  for  the  gas  and  the  electromag¬ 
netic  fields  are  strongly  coupled  due  to  the  presence  of  the 
Joule  heating  term  in  Eq.  (26)  and  the  electrical  conduc¬ 
tivity  in  Eq.  (9).  This  suggests  to  adopt  a  fully  coupled 
approach  by  casting  Eqs.  (6)  and  (23): 


<9rrU  drF  _ 
dt  +lh  ~ 


(27) 


where  conservative  variable,  flux  and  source  term  vectors 
are  now  U  =  (Ug,  Uem),  F  =  (Fg,  Fem)  and  S  = 
(Sg,  Sem).  The  matrix  T  in  Eq.  (27)  reads: 
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(28) 


i  £  if,  j  £  If,  s,p  £  <Sh,  where  5SiPj  =  5sp  Sij,  with  <5 
being  Kronecker’s  delta. 


where  the  translational-rotational  thermal  conductivity  of 

the  heavy-particle  gas  is  Atr  =  At  +  Ar.  3.2.  Spatial  discretization 


Governing  equations 


The  governing  equations  for  the  gas  chemical  composi¬ 
tion  and  temperature  distribution  in  the  ICP  torch  are: 


dr Ug 
dt 
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drF g 
dr 


=  r  Sg 


(23) 


The  application  of  the  Finite  Volume  method  to  Eq.  (27) 
leads  to  the  following  ODE  describing  the  time-evolution 
of  the  conservative  variables  of  cell  c: 

<9U 

r ~gfrcArc  +  rc+\Fc+ 1  -  rc_iFc_i  =  ScrcArc 

(29) 

with  the  cell  volume  (length)  and  centroid  location  be¬ 
ing  Arc  =  rc+ 1/2  -  rc_i/2  and  rc  =  l/2(rc+1/2  + 


rc- 1/2  )>  respectively.  The  evaluation  of  the  diffusive 
flux  Fc+1/2  is  performed  by  approximating  the  values 
and  the  gradients  of  a  given  quantity  p  (e.g.,  tempera¬ 
tures  and  electric  field  components)  by  means  of  an  arith¬ 
metic  average  and  a  second  order  central  finite  differ¬ 
ence,  respectively.  To  facilitate  the  implementation  of 
the  constant  pressure  constraint,  the  solution  update  is 
performed  on  primitive  variables  consisting  of  mass  frac¬ 
tions,  temperatures  and  electric  field  components,  P  = 
(ys,ySi,  T,  Tv,Te ,  Ele,  Eim): 


dP 

rTc-^rcArc  +  rc+ 1  Fc+i  -  rc_  iFc_i  =  ScrcArc 

(30) 

The  transformation  matrix  T  can  be  easily  obtained 
from  the  time-derivative  of  the  conservative  variables 
(dXJ/df)  by  exploiting  the  global  continuity  equation 
(i.e.,  dp/dt  =  0). 


3.3.  Temporal  discretization 


Equation  (30)  is  integrated  in  time  by  means  of  the  back¬ 
ward  Euler  method  [66] : 


4.  COMPUTATIONAL  RESULTS 


The  gas  contained  in  torch  consists  of  nitrogen  molecules 
and  the  related  dissociation  and  ionization  products  (i.e., 
S  =  {e,  N,  N2,  Nj,  N+}).  Simulations  are  performed 
by  means  of  the  ABBA  State-to-State  (StS)  model  [26-29] 
and  the  multi-temperature  (MT)  model  proposed  by  Park 
[43], 

The  torch  radius,  the  number  of  coils  per  unit-length,  the 
frequency  of  the  primary  current  and  the  wall  tempera¬ 
ture  are  set  to  0.08  m,  50,  0.5  MHz  and  350  K,  respec¬ 
tively.  The  current  intensity  of  the  primary  circuit  is  not 
prescribed.  Instead  the  former  is  found  from  the  solution 
by  imposing  that  the  dissipated  power  (per  unit-length)  in 
the  plasma: 

P  =  2n  f  fij  r  dr,  (34) 

Jo 

is  equal  to  Pq  =  350  000W/m.  In  order  to  match  the 
condition  P  =  Po  at  steady-state,  the  current  intensity  is 
multiplied  by  the  scaling  factor  7  =  y7 Pq/P  after  updat¬ 
ing  the  solution  at  each  time-step  [12,  16].  Three  different 
values  are  adopted  for  the  pressure  in  the  torch:  3000  Pa, 
5000  Pa  and  10  000  Pa. 
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where  <5P”  =  P"+1  —  P".  The  local  time-step  At 
is  computed  based  on  the  von  Neumann  number  (£)  as 
At  =  £/(2 pd)  [67],  where  quantity  pd  stands  for  the 
spectral  radius  of  the  diffusive  flux  Jacobian  dF/dJJ. 


In  order  to  advance  the  solution  from  the  time-level  n  to 
the  time-level  n+1,  Eq.  (31)  is  linearized  around  the 
time-level  n.  To  facilitate  the  linearization  procedure,  the 
flux  F  is  written  as  F  =  AdP/dr  and  the  matrix  A 
is  kept  frozen  during  the  linearization.  The  source  term 
is  linearized  by  evaluating  the  related  Jacobian  analyti¬ 
cally  to  enhance  stability.  The  outcome  of  the  lineariza¬ 
tion  procedure  is  a  block-tridiagonal  algebraic  system  to 
be  solved  at  each  time-step: 


Ml"  5P"_!  +  Mc"  5 P’1  +  Mr"  6Pnc+1  =  R"  (32) 
where  the  right-hand-side  residual  is: 

R?;  =  -  (rc+  iFc+i  -rc_iFc_.j  +ScrcArc  (33) 

The  left  (L),  central  (C)  and  right  (R)  block  matrices 
can  be  found  from  the  linearization  procedure  (see,  for 
instance,  [68]  for  more  details).  The  block-tridiagonal 
system  (32)  is  solved  by  means  of  Thomas’  algorithm 
[67]  and  the  solution  updated  at  the  time-level  n  +  1, 

pn+l  =  pn  +  £pn  Xhis 

process  is  continued  until 

steady-state  is  reached. 


In  all  the  cases,  the  solution  is  initialized  with  a  uni¬ 
form  equilibrium  distribution  at  7500  K  (with  both  elec¬ 
tric  field  components  equal  to  0.1  V/m).  The  starting  of 
the  simulation  is  quite  violent.  Severe  temperature  and 
chemical  composition  gradients  occur  at  the  wall,  due  to 
the  cold  temperature  imposed.  During  this  initial  phase 
(in  particular  when  using  the  StS  model),  the  time-step 
must  be  limited  to  avoid  numerical  instabilities.  For  this 
reason,  during  the  first  few  iterations,  the  von  Neumann 
number  is  kept  within  the  range  1000-10  000.  Once  the 
initial  transient  passed,  its  value  is  increased  interactively 
(up  to  500  000  000)  to  accelerate  convergence  to  steady- 
state.  Such  large  values  are  possible  due  to  the  adoption 
of  a  fully  coupled  implicit  time-integration  method. 

For  the  conditions  adopted  in  the  present  work,  the  vibra¬ 
tional  and  free-electron  temperature  profiles  are  essen¬ 
tially  on  top  of  each  other.  This  is  due  to  the  efficient 
energy  exchange  in  N2-e_  interactions  [53] .  For  this  rea¬ 
son,  in  what  follows,  the  vibrational  and  free-electron 
temperatures  are  indicated  as  a  unique  temperature  by 
means  of  the  notation  Tve. 


4.1.  Assessment  of  NLTE  effects 


Before  investigating  in  detail  NLTE  effects,  an  LTE  simu¬ 
lation  is  performed  and  compared  with  the  NLTE  results 
to  assess  the  extent  of  the  departure  from  equilibrium. 
The  pressure  is  set  to  lOOOOPa  as,  for  this  relatively  high 
value,  LTE  conditions  are  often  assumed  [17]. 


Boundary  conditions  are  implemented  through  ghost 
cells  [67]  and  the  related  folding  in  the  central  block  ma¬ 
trices  Me  is  performed  as  suggested  by  Candler  [59]. 


Figure  2  compares  the  LTE  and  NLTE  translational- 
rotational  temperature,  Joule  heating  and  electric  field 
distributions.  Close  to  the  wall,  the  curvature  of  the  LTE 


(a)  (b) 


Figure  2:  LTE  and  NLTE  translational-rotational  temperature  (a),  Joule  heating  (b)  and  induced  electric  field  magnitude 
(c)  distributions  at  10  000  Pa  (unbroken  line  LTE,  dashed  line  MT,  dotted-dashed  line  StS). 


temperature  distribution  changes  sign.  This  is  a  conse¬ 
quence  of  the  non-monotone  behavior  of  the  equilibrium 
total  thermal  conductivity  of  the  working  gas  (nitrogen). 
The  same  trend  is  also  found  in  the  MT  solution  (though 
less  enhanced),  while  the  one  predicted  by  the  StS  model 
does  not  exhibit  any  sign  change  in  its  curvature.  Both 
NLTE  models  predict  that  the  gas  is  in  thermal  equilib¬ 
rium  close  to  the  axis  (the  free-electron  temperature  is 
not  shown  in  Figure  2a),  though  the  equilibrium  is  dif¬ 
ferent  between  the  MT  and  StS  solutions.  In  both  the 
LTE  and  NLTE  simulations,  the  temperature  is  maxi¬ 
mum  on  the  axis,  due  to  the  absence  of  radiative  losses 
in  the  plasma  [7,  8].  Overall,  the  LTE  solution  predicts 
higher  temperature  values,  with  the  difference  being  max¬ 
imum  on  the  axis.  This  is  a  general  trend  observed  in  all 
the  simulations  performed  in  this  work  (and  also  in  the 
multi-dimensional  results  obtained  by  other  investigators 
in  [  17,  20]).  In  NLTE  conditions,  the  temperature  is  lower 
because  the  plasma  is  heated  over  a  wider  region  com¬ 
pared  to  LTE  conditions.  This  is  confirmed  by  the  Joule 
heating  distribution  shown  in  Figure  2b.  The  results  in 


Figure  2  show  that,  even  when  adopting  relatively  high 
pressure  values,  the  LTE  assumption  can  lead  to  a  severe 
overestimation  of  the  gas  temperature  and,  as  a  related 
consequence,  to  a  prediction  of  the  chemical  composi¬ 
tion.  It  is  worth  to  recall  that  in  the  present  work,  the 
effects  of  a  macroscopic  gas  flow  (which  enhance  non¬ 
equilibrium)  are  neglected. 


4.2.  StS  vs  MT 

After  assessing  the  importance  of  NLTE  effects,  the  pre¬ 
dictions  obtained  by  the  StS  and  MT  models  are  com¬ 
pared  in  Figure  3  in  terms  of  temperature  and  N  mole 
fraction.  For  both  the  StS  and  MT  solutions,  decreas¬ 
ing  the  pressure  has  the  effect  of  enhancing  thermal 
non-equilibrium.  In  the  case  of  the  MT  model,  the 
translational-rotational  temperature  is  maximum  on  the 
axis  and  decreases  monotonically  when  approaching  the 
wall. 


10000 
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Figure  3:  Comparison  between  the  StS  and  MT  solutions  in  terms  of  temperature  (left)  and  N  mole  fraction  (right)  at 
different  pressures:  (a)-(b)  p  =  10  000  Pa,  (c)-(d)  p  =  5000  Pa,  (e)-(f)  p  =  3000  Pa  (in  (a),  (c)  and  (e)  unbroken  line  T 
StS,  dashed  line  Tve  StS,  dotted-dashed  line  T  MT,  dotted  line  Tve  MT;  in  (b),  (d)  and  (f)  unbroken  line  StS,  dashed  line 
MT). 
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Figure  4:  Normalized  population  of  the  electronic  levels  of  N  at  different  pressures:  (a)  p  =  10  000  Pa,  (b)  p  =  5000  Pa, 
(c)  p  =  3000  Pa  (line  with  circles  r  =  0  m  [axis],  line  with  squares  r  =  0.04  m,  line  with  triangles  r  =  0.08  m  [wall]). 


On  the  other  hand  the  free-electron  temperature  in¬ 
creases,  reaches  a  maximum  and  then  decreases  till 
reaching  the  value  determined  from  the  adiabatic  bound¬ 
ary  condition.  The  observed  behavior  is  due  to  the  bal¬ 
ance  between  the  Joule  heating  (which  heats  up  the  elec¬ 
tron  gas)  and  the  energy  loss  in  elastic  and  inelastic  col¬ 
lisions,  and  chemical  reactions.  The  peak  location  of  the 
free-electron  temperature  moves  towards  the  wall  when 
decreasing  the  pressure.  This  is  a  consequence  of  the 
Joule  heating  distribution  (not  shown  in  Figure  3)  which 
becomes  sharper  and  clustered  to  the  wall  at  lower  pres¬ 
sures. 

It  is  worth  to  notice  that,  in  the  StS  simulation,  (i)  the 
translational-rotational  temperature  no  longer  exhibits  a 
monotone  behavior  and  (ii)  the  axis  values  of  both  tem¬ 
peratures  are  systematically  lower  than  those  predicted 
by  the  MT  model.  The  differences  observed  between  the 
StS  and  MT  distribution  have  an  effect,  as  it  should  be 
expected,  on  the  chemical  composition  of  the  gas  (see 
Figures  3b,  3d  and  3f). 


Figure  4  shows  the  normalized  population  of  the  elec¬ 
tronic  levels  of  N  on  the  torch  axis  (circles),  in  the  mid¬ 
point  of  the  torch  (squares)  and  at  the  wall  (triangles)  at 
different  pressures.  The  population  distributions  exhibit 
significant  distortions  from  a  Boltzmann  shape  only  close 
to  the  wall  (where  recombination  occurs  and  the  Joule 
heating  is  close  to  its  maximum  value).  Deviations  from 
a  Boltzmann  distributions  are  more  significant  when  in¬ 
creasing  the  pressure  due  to  higher  recombination  (see 
Figures  3b,  3d  and  3f). 


5.  CONCLUSIONS 


A  tightly  coupled  magneto-hydrodynamic  solver  for  the 
study  of  the  weakly  ionized  plasmas  found  in  RF  dis¬ 
charges  has  been  developed.  A  hierarchy  of  thermophys¬ 
ical  models  have  been  added  to  the  solver  to  model  the 
non-equilibrium  effects  in  atomic  and  molecular  plasmas. 
These  include  LTE,  multi-temperature  and  the  more  so- 


phisticated  State-to-State  models.  The  governing  equa¬ 
tions  for  the  flow  and  electromagnetic  fields  have  been 
written  as  a  system  of  time-dependent  conservation-laws. 
Steady-state  solutions  have  been  obtained  by  means  of  an 
implicit  Finite  Volume  Method. 

Results  obtained  by  using  a  multi-temperature  and  State- 
to-State  models  have  shown  that  the  LTE  assumption 
does  not  hold  and  that  its  use  can  lead  to  a  wrong  predic¬ 
tion  of  the  thermo-chemical  state  of  the  gas.  The  analysis 
of  the  translational-rotational  and  free-electron  tempera¬ 
ture  distributions  indicated  that  non-equilibrium  plays  an 
important  role  close  to  the  walls,  due  to  the  combined  ef¬ 
fects  of  Ohming  heating,  and  chemical  composition  and 
temperature  gradients.  The  accurate  investigation  of  the 
population  of  excited  electronic  states  has  shown  that, 
in  view  of  the  absence  of  a  macroscopic  gas  flow,  non- 
Boltzmann  distributions  are  limited  to  a  narrow  region 
close  to  the  torch  wall. 

Future  work  will  focus  on  (i)  including  the  effects  of 
radiation  and  macroscopic  gas  flow  and  (ii)  using  more 
sophisticated  State-to-State  models  to  better  characterize 
NLTE  effects. 
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